In the document, we will demonstrate the effect of label switching on model fitting, different ways to mitigate or eliminate this effect, goodness of fit test and model seletion methods.

Label switching

Effect of label switching

Label switching usually happens during the fitting of mixture models whose components are overlapping. We use the following example to show the effect of label switching.

library(sppmix)
int_surf <- normmix(ps = c(.3, .3, .4),
                    mus = list(c(0.2, 0.2), c(0.3, 0.5), c(.6, .5)),
                    sigmas = list(.01*diag(2), .01*diag(2), .01*diag(2)),
                    lambda = 100,
                    win = square(1))

The intensity surface of the truth model shows that among the three components, two of them are overlapping.

plot(int_surf, main = "Intensity surface for the true model")

You must enable Javascript to view this page properly.

Generate a point pattern from this model and fit this pattern by using est_mix_damcmc().

pp <- rsppmix(int_surf)
plot(pp, int_surf$mus)

fit <- est_mix_damcmc(pp, m=3, truncate = TRUE)
## 0 points are outside W=[0,1]x[0,1]
## Dataset has 101 points
## Preliminaries done. Starting MCMC...

Get the posterior mean after burnin and plot the estimated intensity surface.

postmean <- get_post(fit)
plot(postmean, main = "Posterior intensity surface")
plotmix_2d(postmean, fit$data)

You must enable Javascript to view this page properly.

The plots of intensity surface shows that all components are in the center of the domain. The model doesn’t fitted well.

Diagnositc for label switching

We can detect label switching by using plot_chains(), i.e. trace plot for the MCMC chain of \(\mu\) and p.

plot_chains(fit)

The plots of chains show sudden changes of the parameter values, which indicate label switching. Also, we can use test_labswitch() to detect label switching.

test_labswitch(fit$genmus)
## 
## Checking for label switching...
## Label switching present. 
## Permute the labels to get a better fit,
##           
## or obtain the average of the surfaces
## 
## [1] TRUE

Remedial methods for label swithcing

After label swithcing is detected, we can use FixLS_da() function to mitigate the effect of label switching.

fit2 <- FixLS_da(fit)
postmean2 <- get_post(fit2)
plot(postmean2, main = "Posterior intensity surface after fixing label switching")
plotmix_2d(postmean2, pattern = fit$data)

You must enable Javascript to view this page properly.

Also, to get the intensity surface, we can use function plot_avgsurf() to plot the average of posterior surfaces, which will not be affected by label switching.

plot_avgsurf(fit)

You must enable Javascript to view this page properly.

Goodness of fit test

Function mc_gof() can be used for goodness of fit test. mc_gof() uses Monte Carlo realizations of pattern patterns to simulated the distribution of test statistic and conducts the T-test.

gofts <- mc_gof(pp, postmean2, alpha = 0.05)
## Dataset has 101 points
## Preliminaries done. Starting MCMC...
## 
##          Monte Carlo Goodness of Fit Test 
## Test Statistic: 0.0768  Critical Value: 0.1787 
## Null Hypothesis: Mixture Model fits the pattern well. 
## Alternative: Mixture Model is not appropriate. 
## p-value: 1  Large p-value, Mixture Model fits the data well.
gofts
## NULL

Model selection

We can use function selectMix() to select best number of components and also fit the model by using the best number of component. Notice, all possibel numbers of components need to be specified.

best_fit <- selectMix(pp, 1:5)
best_fit
## $AIC
## [1] -111.0087 -166.1828 -156.1011 -146.0469 -134.9907
## 
## $BIC
## [1]  -95.31798 -134.80136 -109.02898  -83.28401  -56.53713
## 
## $ICLC
## [1]  -95.31798 -126.38897  -70.92187  -19.92264   27.99558
## 
## $BayesFactor
##                numerator 1  numerator 2  numerator 3  numerator 4
## denominator 1 1.000000e+00 1.044637e+14 3.617411e+14 1.221699e+15
## denominator 2 9.572707e-15 1.000000e+00 3.462842e+00 1.169496e+01
## denominator 3 2.764408e-15 2.887802e-01 1.000000e+00 3.377274e+00
## denominator 4 8.185324e-16 8.550689e-02 2.960968e-01 1.000000e+00
## denominator 5 4.394874e-16 4.591046e-02 1.589807e-01 5.369212e-01
##                numerator 5
## denominator 1 2.275378e+15
## denominator 2 2.178153e+01
## denominator 3 6.290073e+00
## denominator 4 1.862471e+00
## denominator 5 1.000000e+00
## 
## $Marginal
## [1] 4.164567e+27 4.350459e+41 1.506495e+42 5.087847e+42 9.475965e+42
## 
## $LogLikelihood
## [1] 61.50435 95.09140 96.05057 97.02345 97.49537

The function selectMix() returns AIC, BIC and ICLC values for each model. Minimum value indicates the best model.

Another approach for model selection is to use birth-death MCMC (Stephens M.(2000)). To fit a Birth-Death MCMC, use function est_mix_bdmcmc. In this case, the parameter m is not the component of the mixture, but the maximum component in a Birth-death MCMC algorithm.

fitbd <- est_mix_bdmcmc(pp, 5)
## 
## Dataset has 101 points, maximum # of components requested is 5
## Preliminaries done. Starting Birth-Death MCMC
tb <- tabulate(fitbd$numcomp, fitbd$maxnumcomp)
mp <- barplot(tb,names.arg=1:fitbd$maxnumcomp,
                 xlab="Number of Components",ylab="Iterations",
                 main="Distribution of the number of components"
                 ,ylim=c(0,1.2*max(tb)))

The histogram of number of components suggests the model with 3 components is the best. Plot the posterior sufrace of model with 3 component.

modelbest <- get_post(fitbd, num_comp = which.max(tb))
plot(modelbest, main = paste("Instensity surface for", which.max(tb), "component model."))
plotmix_2d(modelbest, fit$data)

You must enable Javascript to view this page properly.

Also, we can plot the average surfaces of birth-death MCMC.

plot(fitbd)

You must enable Javascript to view this page properly.